/* 
 * File:   matrizPotencia.c
 * Author: pablo
 *
 * Created on 18 de agosto de 2011, 07:44 PM
 */

#include "mex.h"

/* The computational routine */
void matrizPotencia(double *realBase, double *imBase, double *realSalida, double *imSalida, mwSize m, mwSize n) {
/* Tomo un vector fila base y genero columnas como las potencias de cada elemento de la base*/
    /*m debe ser el tamaño de base y n el número de potencias. salida es nxm*/
int i;
int j;
for(i=0;i<n;i++){
    double auxR=realBase[i];
    double auxI=imBase[i];
    realSalida[m*i]=1;
    imSalida[m*i]=0;
    for(j=1;j<(m-1);j++){
        realSalida[m*i+j]   =   realSalida[m*i+(j-1)] * auxR  - imSalida[m*i+(j-1)]   * auxI;
        imSalida[m*i+j]     =   imSalida[m*i+(j-1)]   * auxR  + realSalida[m*i+(j-1)] * auxI;
    }
    realSalida[m*i+m-1]   =   realSalida[m*i+m-2] * auxR  - imSalida[m*i+m-2]   * auxI;
    imSalida[m*i+m-1]     =  0;
}
}

/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    /* nrhs es el número de inputs desde matlab, *prhs es el puntero a una estructura que los contiene*/
    /* nlhs es el número de salidas y *plhs es el puntero a la estructura que los contendrá*/
    
    double *inBaseReal;               /* 1xN input matrix (fila) */
    double *inBaseIm;               /* 1xN input matrix (fila) */
    mwSize M;                  
    mwSize ncols;                   /* size of matrix */
    double *outMatrixReal;              /* output matrix */
    double *outMatrixIm;              /* output matrix */
    
     /* check for proper number of arguments */
    if(nrhs!=2) {
        mexErrMsgIdAndTxt("MyToolbox:matrizPotencia:nrhs","Two inputs required.");
    }
    if(nlhs!=1) {
        mexErrMsgIdAndTxt("MyToolbox:matrizPotencia:nlhs","One output required.");
    }
    /* check that number of rows in first input argument is 1 */
    if(mxGetM(prhs[0])!=1) {
        mexErrMsgIdAndTxt("MyToolbox:matrizPotencia:notRowVector","Input must be a row vector.");
    }
    
    /*  check that number of cols in second input argument is 1 */
    if(mxGetN(prhs[1])!=1) {
        mexErrMsgIdAndTxt("MyToolbox:matrizPotencia:notColumnVector","Input must be a column vector.");
    }

    /* inputs */
    M = mxGetM(prhs[1]);
    ncols = mxGetN(prhs[0]);
    inBaseReal = mxGetPr(prhs[0]);
    if (mxIsComplex(prhs[0])){
        inBaseIm = mxGetPi(prhs[0]);
    }else{
        inBaseIm=mxGetPi(mxCreateDoubleMatrix(M,ncols,mxCOMPLEX));
    }

    /* create the output matrix */
    plhs[0] = mxCreateDoubleMatrix(M,ncols,mxCOMPLEX);

    /* get a pointer to the real data in the output matrix */
    outMatrixReal = mxGetPr(plhs[0]);
    outMatrixIm = mxGetPi(plhs[0]);

    /* call the computational routine */
    matrizPotencia(inBaseReal, inBaseIm, outMatrixReal, outMatrixIm, M ,ncols);
}

